home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / ASTRNOMY / AA_51.ZIP / TRNSIT.C < prev    next >
C/C++ Source or Header  |  1993-02-09  |  2KB  |  99 lines
  1. /* Calculate time of transit
  2.  * assuming RA and Dec change uniformly with time
  3.  *
  4.  * -- S. L. Moshier
  5.  */
  6.  
  7. #include "kep.h"
  8. extern double glat;
  9.  
  10. /* Earth radii per au */
  11. #define DISFAC 2.3454780e4
  12.  
  13. /* cosine of 90 degrees 50 minutes: */
  14. #define COSSUN -0.014543897651582657
  15. /* cosine of 90 degrees 34 minutes: */
  16. #define COSZEN -9.8900378587411476e-3
  17.  
  18.  
  19. int trnsit(J, lha, dec)
  20. double J; /* Julian date */
  21. double lha; /* local hour angle, radians */
  22. double dec; /* declination, radians */
  23. {
  24. double x, y, z, N, D, TPI;
  25. double lhay, cosdec, sindec, coslat, sinlat;
  26.  
  27. TPI = 2.0*PI;
  28. /* observer's geodetic latitude, in radians */
  29. x = glat * DTR;
  30. coslat = cos(x);
  31. sinlat = sin(x);
  32.  
  33. cosdec = cos(dec);
  34. sindec = sin(dec);
  35.  
  36. x = (J - floor(J-0.5) - 0.5) * TPI;
  37. /* adjust local hour angle */
  38. y = lha;
  39. while( y < -PI )
  40.     y += TPI;
  41. while( y > PI )
  42.     y -= TPI;
  43. lhay = y;
  44. y =  y/( -dradt/TPI + 1.00273790934);
  45. N = x - y;
  46. printf( "approx local meridian transit" );
  47. hms( N );
  48. printf( "UT  (date + %.5lf)\n", N/TPI );
  49.  
  50. if( (coslat == 0.0) || (cosdec == 0.0) )
  51.     goto norise; 
  52.  
  53. /* The time at which the upper limb of the body meets the
  54.  * horizon depends on the body's angular diameter.
  55.  */
  56. switch( objnum )
  57.     {
  58. /* Sun */
  59.     case 0: N = COSSUN; break;
  60.  
  61. /* Moon, elevation = -34' - semidiameter + parallax
  62.  * semidiameter = 0.272453 * parallax + 0.0799"
  63.  */
  64.     case 3:
  65.         N = 1.0/(DISFAC*obpolar[2]);
  66.         D = asin( N ); /* the parallax */
  67.         N =  - 9.890199094634534e-3
  68.             + (1.0 - 0.2725076)*D
  69.             - 3.874e-7;
  70.         N = sin(N);
  71.         break;
  72.  
  73. /* Other object */
  74.     default: N = COSZEN; break;
  75.     }
  76. y = (N - sinlat*sindec)/(coslat*cosdec);
  77. if( (y < 1.0) && (y > -1.0) )
  78.     {
  79. /* Derivative of y with respect to declination
  80.  * times rate of change of declination:
  81.  */
  82.     z = -ddecdt*(sinlat + COSZEN*sindec);
  83.     z /= TPI*coslat*cosdec*cosdec;
  84. /* Derivative of acos(y): */
  85.     z /= sqrt( 1.0 - y*y);
  86.     y = acos(y);
  87.     printf( "rises approx " );
  88.     D = -dradt/TPI + 1.00273790934;
  89.     N = x - (lhay + y)*(1.0 + z)/D;
  90.     hms(N);
  91.     printf( "UT,  sets approx " );
  92.     N = x - (lhay - y)*(1.0 - z)/D;
  93.     hms(N);
  94.     printf( "UT\n" );
  95.     }
  96. norise:        ;
  97. return(0);
  98. }
  99.